% setup env
clear all;
close all;
clc;

% schem parameter
R = 120;
L = 100e-6;
C = 33e-12;

% Resonance friquecy
f0 = 1 / (2*pi * sqrt(L*C));

%Discrete time
T = 10e-5 / f0;
Tmod = 20 / f0;

% time vector
t = 0:T:Tmod;

% Testing model
E = 0*(t<Tmod/4)+1*(t>=Tmod/4);
%plot(t,E);

% init array
U = nan(1, length(t));
dU = nan(1, length(t));
I = nan(1, length(t));

% initial condition
U(1) = 0;
dU(1) = 0;
I(1) = 0;

for k = 2:length(t)
    U(k) = U(k-1) + dU(k-1)*T;
    I(k) = I(k-1) + T*(E(k) - I(k-1)*R -U(k-1)) / L;
    dU(k) = I(k) / C;
end

% plottig
figure;
plot(t,E,t,U);
grid on
title('Step function');
xlabel('t, s');
ylabel('U, V');

% Evaluate AFC
Fmin = 2/Tmod;
Fmax = 2 * f0;
A = 1;
% Friquency array
f = Fmin:((Fmax-Fmin) / 1e+3):Fmax;
K = nan(1, length(f));

U(1) = 0;
dU(1) = 0;
I(1) = 0;

%figure;
for i = 1:length(f)
    E = A * cos(2*pi*f(i)*t);
    
    for k = 2:length(t)
        U(k) = U(k-1) + dU(k-1)*T;
        I(k) = I(k-1) + T*(E(k) - I(k-1)*R -U(k-1)) / L;
        dU(k) = I(k) / C;
    end
    Us = U(fix(end/2):end);
    K(i) = (max(Us)-min(Us)) / 2;
    %plot(t,U);
    %xlabel('t, s');
    %ylabel('U, V');
end

figure;
plot(f,K);
grid on
title('Amplitude-Friquency Characteristic');
xlabel('f, Hz');
ylabel('K');

% Gauss noise

E = 13 * randn(1, length(t));

U = nan(1, length(t));
dU = nan(1, length(t));
I = nan(1, length(t));

U(1) = 0;
dU(1) = 0;
I(1) = 0;

for k = 2:length(t)
    U(k) = U(k-1) + dU(k-1)*T;
    I(k) = I(k-1) + T*(E(k) - I(k-1)*R -U(k-1)) / L;
    dU(k) = I(k) / C;
end

figure;
plot(t,E,t,U);
grid on
title('Gauss noise');
xlabel('t, s');
ylabel('U, V');
